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1.  SUMMARY 


In  this  project  we  developed  a  technique  of  waveform  modeling,  based  on  waveform 
fitting  by  synthetic  seismograms,  and  demonstrated  its  application  to  determine  the  crust 
and  upper  mantle  velocity  structure  beneath  Africa,  China,  and  Canada.  In  our  method  we 
generate  the  synthetic  seismograms  by  the  reflectivity  method,  and  fit  the  observed 
waveforms  by  global  optimization  using  a  Very  Fast  Simulated  Annealing.  Our  technique 
is  complementary  to  the  receiver  function  method  in  that  it  retains  its  advantages,  uses  a 
different  part  of  the  seismogram,  is  sensitive  to  both  P-  and  S-wave  velocities  directly, 
and  obtains  helpful  constraints  in  model  parameters  in  the  vicinity  of  the  Moho.  The 
method  is  particularly  beneficial  if  the  objective  of  using  the  velocity  model  is  to 
determine  the  location  and  focal  depths  of  small,  regional  seismic  events,  due  to  the 
characteristic  feature  of  the  shear-coupled  PL  phase,  which  this  technique  models 
wherever  available,  that  it  samples  a  broader  area  beneath  the  seismic  station  thereby 
representing  a  broad  regional  average.  The  technique  also  inverts  for  the  P-  and  S-wave 
velocities  independent  of  each  other. 

2.  INTRODUCTION 

Fundamental  to  the  study  of  earthquakes  is  their  accurate  location.  Among  the 
location  parameters,  focal  depth  tends  to  be  the  most  unconstrained.  The  primary  reason 
for  this  is  an  inadequate  knowledge  of  local  and  regional  velocity  structure.  Seismologists 
use  two  broad  categories  of  methods  involving  active  and  passive  sources  to  determine 
the  velocity  structure  of  the  Earth.  Active  source  seismology,  however,  is  cost- 
prohibitive,  and  in  many  cases  not  feasible.  Therefore,  passive  source  seismology  which 
employs  natural  earthquakes  as  the  source  to  generate  seismic  waves  is  mostly  used  to 
determine  the  velocity  structure  of  the  Earth.  Nonetheless,  passive  source  seismology  can 
also  be  less  effective  if  the  region  under  consideration  has  relatively  low  rates  of 
seismicity.  Hence,  rather  than  using  common  methods  that  analyze  local  and  regional 
earthquakes,  it  is  necessary  to  develop  and  employ  techniques  that  use  distant  seismic 
sources  to  determine  structure  beneath  such  regions.  In  that  category,  two  commonly 
used  methods  are  receiver  function  analyses  (e.g.,  Owens  et  al.,  1984;  Ammon  et  al., 
1990),  and  surface  wave  dispersion  studies  (e.g.,  Julia  et  al.,  2000;  Pasyanos  and  Walter, 
2002;  Pasyanos,  2005).  These  methods  are  being  used  extensively  in  earthquake 
seismology  to  determine  the  one-dimensional  shear  wave  velocity  structure  and  depth  of 
the  crust-mantle  boundary  (Moho)  beneath  individual  seismic  stations  using  the 
reverberations  of  P-,  and  sometimes  S-wave  phases  from  earthquakes. 

The  objective  of  this  project  was  to  develop  a  third  technique,  based  on  waveform 
fitting  by  synthetic  seismograms,  and  demonstrate  its  application  to  determine  the  crust 
and  upper  mantle  structure  of  the  Earth.  In  our  technique  we  utilize  the  reflectivity 
method  (Kennett,  1983)  to  compute  the  synthetic  seismograms  for  an  earthquake 
recorded  at  a  particular  seismic  station  and  implement  a  global  optimization  algorithm 
using  a  Very  Fast  Simulated  Annealing  (VFSA)  (Ingber,  1989;  Sen  and  Stoffa,  1995)  to 
invert  for  a  one-dimensional  velocity  structure  beneath  that  station.  Our  technique  is 
complementary  to  the  receiver  function  method  in  that  it  retains  its  advantages,  uses  a 
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different  part  of  the  seismogram,  is  sensitive  to  both  P-  and  S-wave  velocities  directly, 
and  obtains  helpful  constraints  in  model  parameters  in  the  vicinity  of  the  Moho.  The 
method  is  particularly  beneficial  if  the  objective  of  using  the  velocity  model  is  to 
determine  the  location  and  focal  depths  of  small,  regional  seismic  events  (Pulliam  et  al., 
2002).  The  primary  reason  for  this  is  the  characteristic  feature  of  the  shear-coupled  PL 
phase,  which  this  technique  models  wherever  available,  that  it  samples  a  broader  area 
beneath  the  seismic  station  thereby  representing  a  broad  regional  average  (Zhao  et  al., 
1996;  Zhao  and  Frohlich,  1996;  Pulliam  et  al.,  2002).  Moreover,  the  technique  also  has 
the  ability  to  invert  for  the  P-  and  S-wave  velocities  independent  of  each  other. 

3.  METHODOLOGY 

3.1  Teleseismic  Phases  Modeled 

In  order  to  improve  upon  the  sampling  area  beneath  a  seismic  station,  relative  to 
receiver  function  inversions,  when  determining  regional  crustal  velocity  structure,  and  to 
obtain  direct  estimates  of  the  P-  and  S-wave  velocities,  we  modeled  the  S,  Sp,  SsPmP, 
and  shear-coupled  PL  (SPL)  phases  in  the  waveforms  of  selected  earthquakes.  In  Figures 
la  and  b  we  show  typical  paths  of  these  waveforms  through  the  crust  and  upper  mantle 
from  a  teleseismic  source.  Correspondingly  in  Figure  lc  we  show  these  waveforms 
generated  synthetically  using  PREM  as  the  velocity  model,  and  for  a  teleseismic  source 
located  at  a  focal  depth  of  600  km  and  recorded  at  an  epicentral  distance  of  50°. 

The  traditional,  direct  S  phase  (path  A-B  in  Figure  la)  is  the  initial,  relatively 
sharp  and  higher  amplitude,  pulse-like  arrival  (Figure  lc)  that  indicates  the  beginning  of  a 
wave  train  with  generally  longer  periods  and  normal  dispersion  (Pulliam  and  Sen,  2005). 
The  particle  motion  of  this  phase  is  essentially  rectilinear  and  all  the  three  components  of 
motion  are  in  phase  with  each  other. 

Upon  impinging  on  the  crust-mantle  boundary  (Moho)  from  below,  a  portion  of 
the  S-wave  converts  to  P-wave  (path  P-N-B  in  Figure  la)  which  then  travels  through  the 
crust  to  arrive  at  the  seismic  station  as  a  precursor  to  the  S-wave  (Figure  lc).  This  phase 
is  called  Sp,  has  been  used  to  model  the  crust  in  earlier  studies  (e.g.,  Jordan  and  Frazer, 
1975;  Owens  and  Zandt,  1997),  and  observed  globally  (Baag  and  Langston,  1986;  Bock, 
1988;  Bock,  1991;  Bock  and  Kind,  1991;  Owens  and  Zandt,  1997).  This  waveform, 
however,  samples  a  more  localized  region  beneath  the  seismic  station  causing  it  to  be  less 
representative  of  a  broader  region  and  more  similar  to  that  of  the  P-coda  (Pulliam  and 
Sen,  2005). 

The  SsPmP  phase  arrives  at  the  base  of  the  crust  as  a  S-wave,  continues  to  travel 
through  the  crust  upward  as  a  S-wave,  but  then  converts  to  a  P-wave  upon  reflecting  off 
the  surface  of  the  Earth,  and  finally  bounces  off  the  Moho  once,  to  arrive  at  the 
seismographic  station  as  a  P-wave  (Langston,  1996;  Owens  and  Zandt,  1997;  Pulliam  and 
Sen,  2005)  (path  Q-M-N-B  in  Figure  la).  Langston  (1996)  demonstrated  that  this 
waveform  is  observed  both  at  regional  and  teleseismic  distances.  He  also  noted  however, 
that  its  propagation  is  dependent  on  the  source  depth  and  distance  because  it  experiences 
a  distance  cross-over  in  the  travel  times  causing  it  to  arrive  before  the  direct  S-wave  at 
regional  distances,  and  after  the  S-wave  at  teleseismic  distances,  with  equal  or  larger 
amplitude  (Figure  lc),  sometimes  interfering  with  and  thus  distorting  the  direct  S-pulse. 
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The  dispersive  wave  train  that  sometimes  follows  S  (Figure  lc)  was  named  the 
“shear-coupled  PL”  (SPL)  phase  by  Oliver  (1961),  who  noted  that  it  is  analogous  to  the 
PL  wave  train  and  appears  after  S  arrivals  at  regional  and  teleseismic  distances.  The 
particle  motion  of  the  SPL  phase  is  prograde  elliptical  in  nature  and  is  confined  to  the 
vertical  plane  (Oliver,  1961).  Based  on  observed  group  and  phase  velocities,  Oliver 
(1961)  also  theorized  that  the  SPL  phase  primarily  is  generated  as  a  S-wave  and  then 
travels  through  the  Earth’s  mantle  until  it  impinges  upon  the  Moho,  thereafter  traveling  as 
trapped  P-waves  and  leaky  SV-waves  (Figure  lb).  This  hypothesis  was  later  validated 
through  presentation  of  computational  methods  for  synthesizing  SPL  by  Chander  et  al. 
(1968),  Frazer  (1977),  and  Baag  and  Langston  (1985).  SPL  phases  were  also  observed  on 
regional  and  teleseismic  seismograms  of  shallow  and  deep  earthquakes  recorded  in 
Europe,  North  America,  central  Andes,  and  Asia  (Zandt  and  Randall,  1985;  Zhang  and 
Langston,  1996;  Zandt  et  al.,  1996;  Owens  and  Zandt,  1997;  Swenson  et  al.,  1999).  The 
SPL  phase  samples  a  broader  region  beneath  the  seismographic  station  (Figure  lb) 
compared  to  the  Sp  phase,  and  its  generation  and  propagation  is  affected  by  seismic 
velocity  gradients,  Vp/Vs  ratios,  impedance  contrasts  across  the  Moho,  and  thicknesses  of 
the  layers  inside  the  Earth  (Baag  and  Langston,  1985;  Pulliam  et  al.,  2002).  With  the 
exception  of  the  direct  S  phase,  all  the  other  phases  discussed  above  appear  prominently 
only  on  the  vertical  and  radial  component  seismograms  at  teleseismic  distances  (Figure 
lc),  due  to  their  conversion  from  S-  to  P-type  motion  at  some  point  in  their  paths. 


Figure  1:  (a)  Typical  ray  paths  for  S  (path  A-B),  Sp  (path  P-N-B),  and  SsPmP  (path  Q- 
M-N-B)  phases  (Modified  from  Pulliam  and  Sen,  2005)  (b)  Propagation  characteristics 
and  excitation  of  shear-coupled  PI  waves  with  distance  and  corresponding  phase  velocity 
(Vph)-period  (T)  curve:  ocn  and  pN  are  the  P  and  S  wave  velocities  at  the  Moho,  and  ci,  C2, 
and  C3  are  the  corresponding  phase  velocities  as  indicated  (Modified  from  Baag  and 
Langston,  1985). 
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Figure  1:  (c)  Three-component  synthetic  seismograms  generated  by  the  parallelized 
reflectivity  program  used  in  this  study,  using  PREM  as  the  velocity  model,  and  for  a 
teleseismic  source  located  at  an  epicentral  distance  of  50°  and  focal  depth  of  600  km.  The 
onset  of  the  S,  Sp,  SsPmP,  and  SPL  phases  are  indicated  on  the  radial  component 
seismogram  (Modified  from  Pulliam  and  Sen,  2005). 

3.2  Waveform  Modeling 

Our  waveform  modeling  technique  combines  a  novel  implementation  of  the 
reflectivity  method  (Kennett,  1983;  Mallick  and  Frazer,  1987)  with  a  global  optimization 
algorithm  (Sen  and  Stoffa,  1991;  1995).  We  compute  the  combined  response  of  all  layers 
of  a  candidate  one-dimensional  Earth  model  using  the  reflectivity  method.  The 
reflectivity  calculation  involves  computation  of  reflectivity  matrices  for  a  stack  of  layers 
as  a  function  of  ray  parameter  and  angular  frequency,  and  produces  all  the  phases 
possible  for  the  specified  stack  of  layers,  source  depth,  and  epicentral  distance.  The 
computations  of  the  reflectivity  responses  for  different  ray  parameters  and  frequencies 
are  completely  independent  of  each  other.  We  use  this  independence  to  adapt  the 
reflectivity  program  to  parallel  computer  architectures,  thereby  decreasing  the 
computation  time  nearly  linearly  with  the  number  of  processors  used  (Pulliam  and  Sen, 
2004).  We  carry  out  message  passing  therein  using  the  Message  Passage  Interface  (MPI) 
Standard  (Gropp  and  Lusk,  1995).  We  distribute  the  computation  over  the  ray 
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parameters,  and  finally  assemble  the  partial  responses,  and  apply  the  inverse 
transformation  from  ray  parameter  to  offset  (a  plane  wave  transformation)  to  generate 
synthetic  seismograms  at  the  required  azimuths  and  distances. 

Following  the  development  of  the  forward  problem,  we  perform  an  optimization 
procedure  to  determine  for  a  given  source-receiver  pair,  the  model  that  produces 
synthetic  waveforms  which  “best  fit”  the  data.  The  criterion  we  use  to  determine  the  best 
fit  is  the  combined  cross-correlation  between  the  vertical  component  of  the  data  and 
synthetics,  and  radial  component  of  the  same,  in  a  specified  time  window.  In  this 
application,  we  define  the  error  as  the  negative  of  a  correlation  function  (Sen  and  Stoffa, 
1991)  given  by: 

E(m)  =  -2[(dv-sv)/{|dv|  +  |sv|}  +  (dr-sr)/{|dr|  +  |sr|}] 
where  dv,  dr,  and  sv,  sr  represent  the  vertical  and  radial  components  of  the  data  and 
synthetics  respectively,  and  |.|  indicates  the  L2  norm. 

Traditionally,  when  the  forward  problem  is  linear  or  there  exists  a  weak  non¬ 
linearity,  derivative-based  methods  such  as  least  squares  are  used  to  solve  the  inverse 
problem  and  estimate  the  model  and  its  uncertainties  (Tarantola,  1994;  Sen  and  Stoffa, 
1995).  However,  in  the  case  of  a  non-linear  problem  such  as  is  common  in  geophysics, 
these  solution  methods  are  generally  not  very  successful.  Therefore,  in  this  method  we 
employ  a  “global  optimization  algorithm”  which  is  only  weakly  dependent  on  the  choice 
of  the  initial  model.  In  particular,  we  use  a  method  called  “Very  Fast  Simulated 
Annealing  (VFSA)”,  which  is  a  variant  of  Simulated  Annealing  (SA)  aimed  at  making 
the  computations  more  efficient  (Ingber,  1989;  Sen  and  Stoffa,  1995). 

Simulated  Annealing  (SA)  is  widely  used  to  attain  a  global,  rather  than  local, 
minimum  while  solving  geophysical  inverse  problems  (Sen  and  Stoffa,  1991;  1995  and 
references  therein).  The  basic  concepts  of  SA  are  derived  from  statistical  mechanics, 
where  an  analogy  is  drawn  between  the  optimization  problem  and  a  physical  system.  SA 
is  analogous  to  the  natural  process  of  crystal  annealing,  in  which  a  solid  in  a  heat  bath  is 
initially  heated  by  increasing  the  temperature  such  that  all  the  particles  are  randomly 
distributed  in  a  liquid  phase,  which  then  gradually  cools.  The  optimization  process 
involves  simulating  the  evolution  of  the  physical  system  as  it  cools  and  anneals  into  a 
state  of  minimum  energy.  At  each  temperature,  the  solid  is  allowed  to  reach  thermal 
equilibrium  where  the  probability  of  it  being  in  that  state  is  given  by  the  Gibbs  or 
Boltzmann  probability  density  function  (Sen  and  Stoffa,  1995).  Very  Fast  Simulated 
Annealing  (VFSA)  is  a  variant  of  SA,  developed  in  order  to  make  it  computationally 
more  efficient.  In  particular,  its  salient  features  include  the  requirement  of  a  temperature 
for  each  model  parameter  which  can  be  different  for  different  model  parameters,  and  the 
use  of  a  temperature  in  the  acceptance  criterion  which  may  be  different  from  the  model 
parameter  temperatures  (Sen  and  Stoffa,  1995).  To  further  illustrate  the  VFSA  technique, 
we  show  a  simplified  flow-chart  in  Figure  2.  The  method  starts  with  an  initial  model  (m°) 
with  an  associated  error  or  energy,  E(m°).  It  then  draws  a  new  model,  mnew,  among  a 
distribution  of  models  from  a  temperature  (T)  dependent  Cauchy-like  distribution,  r(T), 
centered  on  the  current  model  (Figure  2).  The  associated  error  or  energy,  E(mnew),  is  then 
computed  and  compared  with  E(m°)  (Figure  2).  If  the  change  in  energy  (5E)  is  less  than 
or  equal  to  zero,  the  new  model  is  accepted  and  replaces  the  initial  model.  However,  if 
the  above  condition  is  not  satisfied,  mnew  is  accepted  with  a  probability  of  [e5E/T]  (Figure 
2).  This  rule  of  probabilistic  acceptance  in  SA  allows  it  to  escape  a  local  minimum.  We 
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repeat  the  processes  of  model  generation  and  acceptance  a  large  number  of  times  with  the 
annealing  temperature  gradually  decreasing  according  to  a  pre-defined  cooling  schedule 
(Figure  2).  VFSA  is  more  efficient  than  the  traditional  SA  because  it  allows  for  larger 
sampling  of  the  model  space  during  the  early  stages  of  the  waveform  fitting,  and  much 
narrower  sampling  in  the  model  space  as  the  procedure  converges  and  the  temperature 
decreases,  while  still  allowing  the  search  to  escape  from  the  local  minima.  Additionally, 
the  ability  to  perform  different  perturbations  for  different  model  parameters  allows  for 
individual  control  of  each  parameter  and  the  incorporation  of  a  priori  information  (Sen 
and  Stoffa,  1995). 


Figure  2:  Flow-chart  elaborating  the  Very  Fast  Simulated  Annealing  (VFSA)  algorithm 
used  in  this  study  for  the  waveform  inversion  by  global  optimization  (Modified  from  Sen 
and  Stoffa,  1995).  E(m°)  -  error  function  for  the  initial  model  m°.  E(mnew)  -  error 
function  for  the  new  model  mnew,  T  =  temperature,  r(T)  -  temperature  dependent  Cauchy- 
like  distribution. 

3.3  Estimation  of  Uncertainties 

It  is  widely  known  that  solutions  to  geophysical  inverse  problems  are  often  non¬ 
unique.  That  is,  their  error  functions  either  have  broad  minima  or  are  multi-valued, 
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indicating  that  models  that  are  slightly  different  from  the  best-fitting  model  satisfy  the 
data  nearly  as  well,  in  the  first  case,  or  that  one  or  more  very  different  models  also  satisfy 
the  data,  in  the  second  case.  It  is  therefore  necessary  to  explore  the  model  space  and  thus 
identify  the  range  of  models  that  fit  the  data,  and  perhaps  to  identify  characteristics  of  the 
models  that  are  required  by  the  data,  rather  than  which  simply  are  allowed  by  the  data. 
VFSA  conducts  such  a  search  efficiently,  and  the  products  of  multiple  such  searches 
enable  us  to  evaluate  the  uncertainty  in  a  single,  best-fitting  solution.  This  evaluation  is 
particularly  necessary  in  seismic  waveform  modeling  because  more  than  one  model  can 
often  explain  the  observed  data  equally  well  and  trade-offs  between  different  model 
parameters  are  common  (Pulliam  and  Sen,  2005).  The  waveform  modeling  method  we 
developed  incorporates  important  statistical  tools  that  allow  the  user  to  evaluate  the 
uniqueness,  and  physical  feasibility  of  the  resulting  model.  The  most  useful  of  these  tools 
in  evaluating  the  results’  reliability  are  the  Posterior  Probability  Density  (PPD)  function, 
and  the  parameter  correlation  matrix.  To  estimate  these  statistical  parameters  we  cast  the 
inverse  problem  in  a  Bayesian  framework  (e.g.,  Tarantola,  1994;  Sen  and  Stoffa,  1995), 
and  employ  “importance  sampling”  based  on  a  Gibbs’  sampler  (GS)  (Sen  and  Stoffa, 
1995;  Pulliam  and  Sen,  2005).  The  goal  of  “importance  sampling”  is  to  concentrate 
sample  points  in  the  regions  that  are  the  most  “significant”,  in  some  sense  (perhaps,  for 
example,  where  the  error  function  is  rapidly  varying,  or  many  acceptable  solutions  lie). 
Because  this  concentration  is  achieved  using  a  Gibbs’  probability  distribution,  it  has  been 
named  the  “Gibbs’  sampler”  (Sen  and  Stoffa,  1995).  The  Posterior  Probability  Density 
(PPD)  function  [a(m|d„/,.v)]  is  defined  as  a  product  of  a  likelihood  function  [e"£(m)],  and 
prior  probability  density  function,  /?(m).  The  prior  probability  density  function  /Xm), 
describes  the  available  information  on  the  model  without  the  knowledge  of  the  data  and 
defines  the  probability  of  the  model  m  independent  of  the  data.  In  our  applications,  we 
use  a  uniform  prior  within  a  minimum  and  maximum  bound  for  each  model  parameter. 
The  likelihood  function  defines  the  data  misfit  and  its  choice  depends  on  the  distribution 
of  error  in  the  data  (Sen  and  Stoffa,  1995  and  references  therein).  Sen  and  Stoffa  (1996) 
examined  several  different  approaches  to  sampling  models  from  the  PPD  and  concluded 
that  a  multiple- VFSA  based  approach,  though  theoretically  approximate,  is  the  most 
efficient.  In  a  multiple-VFSA  approach  we  make  several  VFSA  runs  (20  in  this  study) 
with  different  random  starting  models  and  use  all  the  models  sampled  along  to 
characterize  uncertainty  in  the  model.  We  use  all  these  sampled  models  to  compute 
approximate  marginal  PPD  and  posterior  correlation  matrices  to  characterize 
uncertainties  in  the  derived  results.  The  posterior  correlation  matrix  measures  the  relative 
trade-off  between  individual  model  parameters  and  is  computed  by  normalizing  the 
covariance  between  two  model  parameters  (Sen  and  Stoffa,  1996).  Computationally,  the 
correlation  between  ith  and  jth  model  parameters  is  given  by  their  covariances  divided  by 
the  square  root  of  the  product  of  the  covariances  of  each  parameter  with  itself.  In  a  later 
section  during  discussion  of  the  application  of  the  technique  to  seismological  data 
recorded  in  Africa,  we  provide  descriptions  of  interpretations  of  the  resulting 
computations  of  the  PPD  and  correlation  matrix. 
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4.  RESULTS  -  APPLICATION  OF  THE  METHOD 


We  applied  our  modeling  method  described  above  to  data  from  large -magnitude, 
deep-focus  earthquakes  recorded  teleseismically  during  1976  -  2005  at  permanent 
broadband  seismic  stations  spanning  the  continent  of  Africa,  China  and  Canada.  The 
focal  depths  of  these  earthquakes  range  between  200  km  and  600  km,  and  their 
magnitudes  lie  between  5.5  and  7.0.  Since  the  goal  is  to  also  model  the  SPL  phase  that  is 
generated  close  to  the  seismic  station  (within  an  area  of  -100  km  x  100  km)  (Frazer, 
1977),  we  chose  such  a  focal  depth  range  to  eliminate  the  SPL  phase  generated  at  the 
earthquake  source.  Epicentral  distances  from  the  seismic  stations  of  the  selected 
earthquakes  are  between  30°  and  80°  so  as  to  avoid  possible  incorporation  of  phases  that 
interacted  with  the  Earth’s  core.  Initially,  we  filtered  the  raw  data  obtained  from  the 
global  database  for  the  selected  events,  using  a  six-pole  Butterworth  bandpass  filter  with 
comer  frequencies  of  0.005  and  0.25  Hertz  respectively.  We  then  decimated  the  data  such 
that  the  sample  interval  is  0.5  seconds.  The  data  window  we  analyzed  includes  30 
seconds  prior  to  the  arrival  of  the  direct  S  phase  and  180  seconds  following  it.  The  choice 
of  the  beginning  time  for  our  data  window  follows  from  a  study  by  Jordan  and  Frazer 
(1975)  who  showed  that  for  a  deep  focus  event  (-600  km)  of  intermediate  magnitude 
(-6.1),  at  teleseismic  distances,  the  Sp  phase  resulting  from  a  single  conversion  at  the 
Moho  (-35  km-40  km),  precedes  the  S  phase  by  about  5-6  seconds.  Since  the  events  we 
modeled  in  these  applications  also  lie  in  that  category,  and  the  only  phase  arriving  before 
the  S  phase  that  we  model  is  the  Sp,  we  do  not  expect  to  observe  any  Moho-converted  Sp 
phase  before  -15  seconds  from  the  S  phase.  Therefore,  the  start  time  of  our  data  window 
(30  seconds  before  the  S  arrival)  provides  ample  lead  time  for  us  to  never  miss  the 
observation  of  the  Sp  phase  if  any.  For  each  station,  the  initial  model  we  chose  is  one 
obtained  from  a  previous  published  study,  wherever  available,  or  Preliminary  Reference 
Earth  Model  (PREM)  (Dziewonski  and  Anderson,  1981)  where  an  earlier  study  has  not 
been  published.  We  also  experimented  with  an  initial  model  consisting  of  crustal  layers 
of  equal  thickness  and  increasing  velocities,  superimposed  on  PREM.  However,  the  final 
models  obtained  using  our  method  were  similar  within  one  standard  error,  thereby 
emphasizing  minimal  dependence  of  our  method  on  the  starting  model.  For  use  in 
reflectivity  computations,  we  also  incorporate  the  source  mechanism  of  each  earthquake 
from  the  Global  CMT  catalog,  and  use  a  Gaussian  source-time  function.  Following 
similar  forward  calculations  for  each  source  -  receiver  pair,  we  carried  out  the  waveform 
fitting  procedure  for  each  using  200  iterations.  Prior  to  our  choice  of  the  number  of 
iterations,  we  experimented  with  200,  400,  600,  and  800  iterations,  and  have  consistently 
observed  that  after  -165  iterations  the  error  reaches  an  optimal  value  and  does  not  change 
with  subsequent  iterations.  This  feature  is  a  diagnosis  in  our  method  to  confirm  that  the 
process  has  converged.  Therefore,  we  chose  200  iterations  as  a  threshold  for  all  our 
computations.  Additionally,  earlier  studies  (Sen  and  Stoffa,  1995  and  references  therein) 
have  shown  that  Very  Fast  Simulated  Annealing  (VFSA)  typically  converges 
significantly  faster  than  other  methods  in  the  category,  hence  the  name.  Based  on 
examples  documented  by  Sen  and  Stoffa  (1995),  we  chose  an  initial  temperature  of  10~3 
units  at  the  start  of  our  waveform  fitting  process  for  each  model  parameter  and  allowed  it 
to  cool  down  to  10~10  units  throughout  the  process.  In  our  method,  we  allowed  each 
model  parameter  (velocity  of  the  P-wave,  Vp,  velocity  of  the  S-wave,  Vs,  thickness  of  the 
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layer,  and  density)  to  vary  within  ±10%  of  initial  values.  We  conducted  trial  runs  with  the 
model  parameters  varied  within  ±10%,  ±15%,  ±20%,  and  ±30%.  Our  results  produced 
similar  final  models  that  were  within  one  standard  error.  Additionally,  a  significant 
variation  in  model  parameters  is  not  realistic  given  the  tectonic  and  geologic  setting  of 
the  regions.  Therefore,  to  maintain  reasonable  computational  time  and  to  allow  variations 
that  are  more  realistic,  we  varied  the  model  parameters  ±10%. 

Below,  we  report  results  of  waveform  fitting  for  selected  teleseismic  earthquakes 
recorded  in  Africa,  China  and  Canada.  Work  is  currently  ongoing  on  the  data  from  China 
and  Canada.  The  detailed  study  of  crust  and  upper  mantle  velocity  structure  beneath 
Africa  using  our  method  is  presented  in  Gangopadhyay  et  al.  (2007),  therefore  herein  we 
only  discuss  important  results  from  that  study.  For  the  seismic  stations  that  recorded 
better  quality  data  overall,  we  show  the  waveform  correlations  for  events  recorded  at  that 
station,  and  also  describe  the  interpretations  of  the  uncertainty  computations  as  an 
example.  A  comment  on  amplitude  matches:  The  most  successful  match  between 
synthetics  and  data  would  be  one  in  which  the  synthetic  waveform  matched  the  data 
exactly  -  wiggle  for  wiggle.  This  is  unrealistic  for  several  reasons,  including  the  fact  that 
models  used  to  compute  synthetics  are  layered,  isotropic,  limited  to  ten  to  sixteen  layers, 
and  have  fixed  attenuation  (Q)  values.  Further,  the  source  time  function  is  assumed  to  be 
Gaussian  and  its  focal  mechanism  is  assumed  to  be  correctly  represented  by  Harvard’s 
CMT  solution.  To  minimize  complexities  in  the  source  time  function  we  avoid  very  large 
earthquakes.  Given  the  uncertainties  in  model  Q  and  focal  mechanisms,  which  will 
largely  control  relative  amplitudes  of  various  phases,  we  focus  our  fitting  criteria  on 
matching  each  phase’s  arrival  time,  polarity,  and  pulse  character.  Fitting  the  amplitude  of 
each  phase,  while  desirable,  is  deemed  to  be  of  lesser  importance. 

4.1  Africa 

We  applied  our  modeling  method  to  data  recorded  at  twelve  permanent  broadband 
seismic  stations  spanning  the  continent  of  Africa  (Figure  3).  A  total  of  seventeen 
earthquakes  were  used  in  this  study  selected  from  the  global  catalog  of  Centroid  Moment 
Tensors  (CMT)  (1976  -  2004).  The  region  encompassing  north  and  west  Africa  includes 
the  seismic  stations  of  TAM,  DBIC,  MBO,  and  MDT  (Figure  3).  Among  these  stations, 
TAM  recorded  data  of  better  quality;  examples  of  the  waveform  fits  for  three  events 
recorded  at  TAM  are  shown  in  Figure  4a.  We  observe  S,  Sp,  and  SsPmP  phases 
consistently  on  all  these  event  seismograms  and  they  correlate  well  with  the  synthetics 
generated  by  the  optimization  technique  (Figure  4a).  On  event  3,  we  also  observe  a 
prominent  SPL  phase  and  the  synthetics  match  it  well  (Figure  4a).  Particle  motion 
diagrams  for  the  corresponding  time  window  on  both  data  and  synthetics  confirm  this 
observation  (Figure  4b).  We  observe  prograde  elliptical  particle  motion,  which  is 
diagnostic  of  the  SPL  phase,  on  both  diagrams  (Figure  4b).  Except  for  event  3,  none  of 
the  others  have  any  signature  of  the  SPL  phase,  as  confirmed  by  particle  motion  diagrams 
for  corresponding  time  windows.  This  observation  -  that  one  source-receiver  pair  would 
show  SPL  while  other,  similar  paths  would  not  -  is  unexpected  and  we  cannot  explain  it. 

Based  on  the  waveform-fitting  results  for  all  six  events  recorded  at  TAM,  we 
generate  P-  and  S-wave  velocity  models  up  to  a  depth  of  100  km  (Figure  4c).  We  observe 
some  variability  in  the  models  (Figure  4c),  and  hence  compute  the  uncertainties  for  each 
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Figure  3:  Map  showing  locations  of  the  earthquakes  used  in  this  study  (blue  stars)  and 
the  permanent  broadband  seismic  stations  in  the  African  continent  (red  triangles)  that 
have  recorded  these  earthquakes.  The  respective  station  codes  are  shown  adjacent  to  each 
seismic  station. 
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Figure  4:  (a)  Vertical  and  radial  component  seismograms  for  three  events  recorded  at 
TAM  showing  the  observed  (solid  line)  and  synthetic  (dashed  line)  waveforms.  The 
correlated  waveforms  are  indicated  on  the  panels,  (b)  Particle  motion  diagrams  for  a  10-s 
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time  window  around  the  SPL  arrivals  for  event  3  showing  the  diagnostic  prograde 
elliptical  motion  of  the  SPL  phase.  The  dotted  portions  of  the  diagrams  indicate  the 
beginning  of  the  motion. 
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Figure  4:  (c)  P-  and  S-wave  velocity  models  up  to  100  km  for  station  TAM  from  the 
inversion  results  for  individual  events  recorded  at  TAM.  (d)  Model  parameter  correlation 
matrices  for  events  1  (left-hand  panel)  and  3  (right-hand  panel)  recorded  at  TAM.  Each 
small  square  represents  a  model  parameter  (Vp,  Vs,  Thickness  of  layer,  and  Density)  on 
both  the  horizontal  and  vertical  axes.  The  correlations  range  between  -1  and  1.  Sparse 
colored  squares  off-diagonal  in  the  lower  crust-upper  mantle  in  event  3  (right-hand  panel) 
compared  to  that  in  event  1  (left-hand  panel)  indicate  better  resolution  and  confidence 
(less  trade-off)  in  this  region. 

model  using  the  statistical  tools  described  earlier  to  choose  the  “best”  model.  In  Figure 
4d,  we  show  examples  of  parameter  correlation  matrices  computed  from  the  modeling 
results  of  events  1  and  3.  Each  small  square  along  an  axis  of  the  parameter  correlation 
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matrix,  either  horizontally  or  vertically,  represents  a  model  parameter  (Figure  4d).  Since 
every  model  layer  consists  of  four  independent  model  parameters  (Vp,  Vs,  thickness,  and 
density),  four  small  squares  combined  together  represent  a  model  layer  on  both  axes 
(Figure  4d).  Correlation  values  range  between  -1  and  1  and  are  symmetric  about  the 
diagonal  of  the  matrix,  hence,  for  clarity,  we  show  only  values  below  the  diagonal 
(Figure  4d).  Values  along  the  diagonal  are  ones,  simply  indicating  that  each  parameter  is 
perfectly  correlated  with  itself.  Off-diagonal  colored  squares  indicate  significant  cross¬ 
correlation  (trade-offs)  between  corresponding  model  parameters.  In  the  parameter 
correlation  matrices  for  both  events  (Figure  4d),  layers  comprising  the  upper  crust  have 
greater  independence,  as  indicated  by  the  sparse  distribution  of  off-diagonal  cross¬ 
correlations  whose  absolute  values  are  greater  than  ±0.5  (colored  squares).  Also,  for  both 
events,  the  level  of  tradeoffs  among  model  parameters  in  these  shallow  layers  is  similar 
(Figure  4d).  For  event  1,  however,  the  layers  comprising  the  lower  crust  and  upper 
mantle  have  larger  off-diagonal  cross-correlations,  indicating  significant  tradeoffs  (Figure 
4d).  On  the  contrary,  for  event  3,  even  the  lower  crustal  and  upper  mantle  layers  appear 
better  constrained  (Figure  4d).  Intriguingly,  the  SPL  phase  is  also  observed  in  the 
seismogram  of  event  3  but  not  in  that  of  event  1  (Figure  4a).  This  observation  attests  to 
the  fact  that,  if  SPL  is  present  in  the  seismogram  and  is  well  modeled,  we  are  able  to 
better  constrain  the  structure  of  the  lower  crust  and  upper  mantle.  This  result,  which  is 
expected,  due  to  the  sensitivity  of  SPL  to  those  parts  of  the  model  (Figure  lb),  drives  our 
decision  to  generate  velocity  models  down  to  the  Moho  for  seismic  stations  at  which  SPL 
is  not  observed. 

We  observed  and  successfully  fit  S  and  SP  phases  in  seismograms  recorded  at 
DBIC  and  at  MBO.  The  SsPmP  phase  appears  on  the  vertical-component  seismograms  of 
an  event  at  DBIC  and  one  at  MBO  but  is  absent  on  their  radial-component  seismograms. 
We  do  not  observe  the  SPL  phase  in  seismograms  recorded  at  DBIC  but  do  so  in  the 
seismograms  at  MBO.  The  SPL  phase  is  prominent  on  both  vertical  and  radial  component 
seismograms  of  the  event  recorded  at  MBO,  which  we  confirm  with  particle  motion 
diagrams  that  show  prograde  elliptical  motion  for  the  corresponding  time  window.  At 
MDT,  within  the  time -window  expected  to  contain  the  phases  analyzed  in  this  study,  we 
were  unable  to  clearly  identify  them  and  they  appear  to  be  contaminated  by  interfering 
arrivals,  and  hence  are  also  not  well-correlated  with  the  synthetics  generated  by  the 
waveform  modeling  process.  Station  MDT  is  located  near  the  Atlas  Mountains  in 
Morocco,  and  so  may  be  underlain  by  complicated  three-dimensional  structure,  which  the 
waveform  modeling  program  used  in  this  study  is  unable  to  model  accurately.  The  lack  of 
proper  correlation  between  the  synthetic  seismograms  and  the  data  at  MDT  is  likely  a 
consequence  of  this  limitation. 

East  Africa  contains  the  permanent  broadband  seismic  stations  ATD,  FURI, 
KMBO,  and  MBAR  (Figure  3).  However,  these  stations  are  situated  within  and  on  the 
flanks  of  the  active  East  African  rift  system  and  therefore  waves  recorded  at  these 
stations  sample  the  complicated  three-dimensional,  anisotropic  structure  beneath  the  rift 
(Ayele  et  al.,  2004;  Dugda  and  Nyblade,  2006).  The  three-dimensional  structure  is 
manifested  in  the  seismograms  as  numerous,  possibly  scattered,  refracted,  or  split,  phase 
arrivals  with  strong  interference  amongst  themselves.  Anisotropy  inferred  from  shear- 
wave  splitting  studies  has  also  been  reported  for  stations  ATD,  FURI,  and  KMBO  by 
Ayele  et  al.  (2004).  As  we  note  earlier,  the  waveform  inversion  method  we  use  in  this 
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study  is  capable  of  estimating  azimuthally  dependent  one-dimensional  (although 
azimuthally  dependent)  structure  only,  hence  waveforms  for  events  recorded  at  these 
stations  are  not  precisely  correlated  in  some  cases.  On  both  the  vertical  and  radial 
component  seismograms  of  most  events  at  ATD,  FURI,  KMBO,  and  MBAR,  we  observe 
S  and  SsPmP  phases.  On  the  other  hand,  except  at  MBAR,  we  observe  the  SP  phase  only 
on  the  vertical  component  for  these  events.  Only  at  FURI,  for  event  12,  and  MBAR,  for 
event  15,  did  we  see  SPL  phase  in  the  seismograms. 

Seismic  stations  TSUM,  LSZ,  LBTB,  and  SUR  are  located  in  southern  Africa 
(Figure  3).  However,  for  the  lone  event  recorded  at  SUR,  we  are  able  to  identify  only  the 
direct  S  phase,  and  thus  we  do  not  attempt  to  generate  P-  and  S-wave  velocity  models  for 
SUR.  We  observe  and  successfully  fit  S,  Sp,  and  SsPmP  phases  on  events  13  and  14 
recorded  at  TSUM,  event  9  recorded  at  LSZ,  and  event  1 1  at  LBTB.  But,  except  for  event 
13  at  TSUM,  we  are  able  to  identify  the  Sp  phase  on  both  the  vertical  and  radial 
component  seismograms  at  TSUM,  LSZ,  and  LBTB.  Similarly,  the  SsPmP  phase  is  only 
identifiable  on  the  vertical  component  seismogram  for  events  14  and  9  recorded  at  TSUM 
and  LSZ  respectively.  It  is,  however,  absent  on  the  seismograms  for  event  13  at  TSUM. 
We  do  not  observe  the  SPL  phase  in  the  seismograms  of  any  event  recorded  at  these 
stations. 

For  each  seismic  station  in  the  three  broad  regions  of  the  African  continent,  based 
on  the  waveform  correlations  described  earlier,  we  generate  azimuthally  dependent  P- 
and  S-wave  velocity  models  beneath  the  station.  At  the  seismic  stations  where  the  SPL 
phase  is  not  observed  and  modeled,  we  only  generate  velocity  models  up  to  the  crust- 
mantle  boundary  (Moho),  since  the  body-wave  phases  Sp,  SsPmP,  and  S  do  not  constrain 
models  below  the  Moho.  Where  SPL  is  observed  we  generate  the  velocity  models  up  to 
an  arbitrarily  chosen  depth  of  100  km.  Wherever  available,  we  compare  our  velocity 
models  with  those  obtained  from  earlier  studies  for  the  stations. 

Figures  5a-d  shows  the  obtained  P-  and  S-wave  velocity  models  for  the  north  and 
west  African  stations  of  TAM,  DBIC,  MBO,  and  MDT.  The  estimates  of  crustal 
thickness  beneath  these  stations  range  between  36  km  and  42  km  (Figures  5a-d), 
comparable  to  regional  estimates  of  34  km  to  40  km  by  Pasyanos  et  al.  (2004).  We  also 
note  that  the  crust  is  slightly  more  thick  in  west  Africa  beneath  seismic  stations  DBIC 
and  MBO  (-41^12  km)  (Figures  5b  and  5c),  compared  to  the  seismic  stations  TAM  and 
MDT  in  north  Africa  (~36  km-38  km)  (Figures  5a  and  5d).  A  similar  observation  was 
also  made  earlier  by  Pasyanos  and  Walter  (2002)  using  surface  wave  dispersion 
tomography,  and  by  Marone  et  al.  (2003)  using  joint  inversion  of  local,  regional,  and 
teleseismic  data.  Except  for  TAM,  the  crust  below  all  the  stations  appears  to  be  fairly 
simple  in  structure  (Figures  4a-d)  (Sandvol  et  al.,  1998),  suggesting  that  it  is  minimally 
affected  by  large-scale  tectonic  processes.  However,  a  middle  to  lower  crustal  low- 
velocity  zone  obtained  beneath  all  the  seismic  stations  in  the  region  (Figures  5a-d), 
indicate  possible  local  tectonic  influences.  Our  estimate  of  crustal  thickness  beneath 
TAM  (~36  km)  is  similar  to  that  obtained  from  receiver  function  studies  by  Sandvol  et  al. 
(1998)  (38  ±  2  km)  (Figure  5a),  from  Rayleigh  wave  group  velocity  dispersion  studies  by 
Hazier  et  al.  (2001)  (43  ±  5  km),  and  from  surface-wave  dispersion  tomography  by 
Pasyanos  and  Walter  (2002)  (~40  km).  TAM  is  close  to  the  location  of  the  Hoggar  hot 
spot  but,  as  noted  by  Sandvol  et  al.  (1998),  the  crustal  thickness  indicates  that  a  mantle 
plume  has  not  significantly  altered  the  crust  here.  However,  in  contrast  to  the  model  of 
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Figure  5:  Preferred  P-  and  S-wave  velocity  models  for  (a)  TAM  (b)  DBIC  (c)  MBO  and 
(d)  MDT.  The  P-  and  S-wave  velocity  models  (broken  lines)  in  (a),  (b),  and  (d)  are  from 
receiver  function  studies  by  Sandvol  et  al.  (1998).  The  dotted  lines  in  (a)  show  the  P-  and 
S-wave  velocity  model  for  TAM  within  ±  two  standard  errors. 

Sandvol  et  al.  (1998),  the  crustal  P-  and  S-wave  velocities  we  obtain  in  this  study  at 
TAM,  are  both  slightly  lower  (Figure  5a).  These  velocities  at  TAM  range  between  6.25 
km/s-6.8  km/s,  and  3.1  km/s-3.9  km/s,  respectively  (Figure  5a).  Upper  mantle  P-wave 
velocities  exhibit  a  gradational  increase  with  depth  below  the  Moho,  whereas  the  S-wave 
velocities  are  nearly  constant  (4  km/s-4.2  km/s)  within  that  range  of  depths  (Figure  5a). 
Furthermore,  we  also  obtain  an  anomalous  low-velocity  zone  of  ~5  km  thickness  at  the 
base  of  the  upper  crust  (Figure  5a)  that  appears  to  be  well  constrained  based  on  the 
uncertainty  estimates  described  earlier  (Figure  4d).  However,  we  are  unable  to  confirm 
the  existence  of  this  layer  from  any  other  independent  studies.  At  the  west  African  coastal 
station  DBIC,  we  obtain  a  Moho  depth  of  ~41  km  which  is  similar  to  that  obtained  by 
Sandvol  et  al.  (1998)  (—40  ±  2.3  km)  (Figure  5b).  P-wave  velocities  range  between  6.7 
km/s-7  km/s  in  the  upper  crust,  and  6.5  km/s-7.3  km/s  in  the  lower  crust  (Figure  5b).  On 
the  contrary,  S-wave  velocities  show  a  gradational  increase  in  the  crust  with  depth  from 
3.7  km/s-4.7  km/s  (Figure  5b).  Anomalous  Vp/Vs  ratios  are  thus  caused  by  a  low  P-wave 
velocity  zone  of  ~15  km  thickness  in  the  lower  crust.  However,  due  to  the  trade-offs 
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between  the  model  parameters  in  this  depth  range,  we  conclude  that  this  anomaly  is  not 
well-constrained.  Beneath  MBO,  any  prior  P-  and/or  S-wave  velocity  models  are  absent. 
Sandvol  et  al.  (1998)  had  analyzed  seismic  data  recorded  at  MBO  but  the  anomalous  data 
did  not  allow  them  to  estimate  a  velocity  model.  Thus,  the  velocity  model  for  MBO  that 
we  obtain  from  this  study,  is  to  our  knowledge,  the  first  estimate  of  its  kind.  We  obtain  a 
crustal  thickness  of  -42  km  (Figure  5c),  which  is  similar  to  that  found  underneath  other 
stations  in  the  region.  A  regional  crustal  thickness  of  43  ±  5  km  obtained  from  Rayleigh 
wave  group  velocity  dispersion  studies  (Hazier  et  al.,  2001)  correlates  well  with  the 
results  of  this  study  at  MBO.  The  P-  and  S-wave  velocities  at  MBO  in  the  crust  range 
between  5.6  km/s-7.2  km/s,  and  3.1  km/s-4.1  km/s,  respectively  (Figure  5c).  We  also 
observe  an  anomalous  lower  crustal,  approximately  15  km-thick  zone  of  relatively  low 
P-  and  S-wave  velocities  (6.8  km/s  and  3.8  km/s)  beneath  MBO  (Figure  5c).  Our  P-  and 
S-wave  velocity  models  beneath  MDT  predict  a  Moho  depth  of  -38  km  (Figure  5d). 
Surprisingly,  even  with  the  poor  waveform  fit  by  synthetics  to  the  event  recorded  at 
MDT,  this  result  is  consistent  with  the  estimates  obtained  by  Sandvol  et  al.  (1998)  (36  ± 

1.3  km).  As  also  noted  by  Sandvol  et  al.  (1998)  and  Pasyanos  and  Walter  (2002),  the 
slightly  shallower  Moho  at  MDT,  compared  to  that  at  DBIC  and  MBO,  indicates  that  in 
spite  of  its  proximity  to  the  Atlas  Mountains,  there  is  no  crustal  thickening  associated 
with  them,  and  a  significant  root  is  absent  beneath  the  mountains.  Sandvol  et  al.  (1998) 
concluded  that  this  may  be  a  possible  outcome  of  the  fact  that  there  existed  a  failed  rift 
earlier  which  was  subsequently  inverted.  In  this  study,  beneath  MDT,  we  obtain  average 
P-  and  S-wave  velocities  in  the  crust  of  -6.4  km/s  and  3.7  km/s,  respectively,  except  in  a 
low-velocity  zone  of  -10  km  thickness  in  the  lower  crust  where  these  are  6.2  km/s  and 

3.4  km/s,  respectively  (Figure  5d).  A  similar  zone  has  also  been  postulated  by  the  earlier 
velocity  model  obtained  from  receiver  function  studies  (Sandvol  et  al.,  1998).  However, 
because  of  poor  waveform  fits  at  MDT  in  this  study,  we  are  unable  to  postulate  the 
existence  of  this  zone. 

For  east  Africa,  we  generate  P-  and  S-wave  velocity  models  beneath  seismic 
stations  ATD,  FURI,  KMBO,  and  MBAR,  which  are  shown  in  Figures  6a-d.  Located 
within  and  on  the  flanks  of  the  East  African  Rift  System,  which  is  relatively  well  studied, 
these  stations  are  situated  in  active  tectonic  domains.  Except  beneath  ATD  (Figure  6a), 
where  the  crust  is  significantly  thin  compared  to  the  other  stations  in  the  region  (Figures 
6d-f),  the  Moho  is  generally  between  -38  km-41  km  deep.  The  estimate  of  crustal 
thickness  beneath  ATD,  however,  is  the  subject  of  an  active  debate.  Using  a  grid  search 
method  to  model  receiver  functions  for  eleven  earthquakes  recorded  at  ATD,  Sandvol  et 
al.  (1998)  obtained  a  crustal  thickness  of  -10  km  (Figure  6a).  But,  Dugda  and  Nyblade 
(2006)  used  H-k  analysis  of  receiver  functions  and  predicted  a  crustal  thickness  of  -23  ± 

1.5  km  beneath  ATD,  consistent  with  earlier  results  from  inversion  of  gravity  data  for  the 
general  area  by  Tiberi  et  al.  (2005).  In  this  study,  our  velocity  model  beneath  ATD  shows 
comparable  velocity  discontinuities  at  -10  km  and  -21  km  depths  (Figure  6a),  suggesting 
that  either  of  these  depths  could  be  interpreted  as  the  Moho.  However,  the  layer  at  10  km 
depth  appears  to  be  poorly  constrained  compared  to  the  layer  at  -21  km  depth,  as 
evidenced  from  the  PPD  (Figure  6b)  and  the  parameter  correlation  matrix  (Figure  6c) 
computations.  Therefore,  we  prefer  a  crustal  thickness  of  -21  km.  Irrespective  of  the 
debate  on  the  crustal  thickness  at  ATD,  the  crust  beneath  it  is  significantly  thinner  than 
the  crust  beneath  other  seismic  stations  in  east  Africa  (Figures  6d-f).  Located  within  the 
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Afar  depression,  close  to  the  coast  of  the  Red  sea  on  the  eastern  edge  of  the  African 
continent,  such  a  thin  crust  is  expected  at  ATD  because  of  highly  stretched  continental 
crust  (Sandvol  et  al.,  1998). 


Figure  6:  (a)  Obtained  velocity  model  from  this  study  (solid  line)  up  to  30  km  for  station 
ATD  and  available  velocity  model  from  Sandvol  et  al.  (1998)  (broken  line).  Posterior 
Probability  Distribution  (b)  and  parameter  correlation  matrix  (c)  for  event  13  recorded  at 
ATD  showing  the  tradeoffs  between  different  model  parameters  in  different  layers,  (d) 
Obtained  velocity  model  from  this  study  (solid  line)  up  to  100  km  for  station  FURI  and 
available  velocity  model  from  Ayele  et  al.  (2004)  (broken  line).  Preferred  P-  and  S-wave 
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velocity  models  for  (e)  KMBO  up  to  a  depth  of  50  km,  and  (f)  MBAR  up  to  a  depth  of 
100  km. 

The  P-  and  S-wave  velocities  in  the  crust  at  ATD  range  between  4.7  km/s-7.2 
km/s,  and  2.5  km/s-4.3  km/s,  respectively  (Figure  6a).  Crustal  P-wave  velocities  below 
about  5  km  depth  are  relatively  high  and,  as  noted  by  Dugda  and  Nyblade  (2006),  could 
indicate  a  highly  mafic  composition  caused  by  igneous  rock  emplacement  during  the  syn- 
rift  stage.  Figure  6d  shows  the  preferred  P-  and  S-wave  velocity  models  beneath  FURI 
which  is  situated  in  the  northern  part  of  the  western  Ethiopian  plateau.  We  obtain  an 
estimate  of  crustal  thickness  beneath  FURI  of  -39  km,  which  is  similar  to  that  obtained 
using  receiver  function  analyses  by  Ayele  et  al.  (2004)  (-40  km),  and  Dugda  et  al.  (2005) 
(-44  km).  The  —40  km  thick  crust  beneath  FURI,  which  is  located  close  to  the  border  of 
the  western  Ethiopian  plateau  and  the  Afar  depression,  is  also  consistent  with  previous 
refraction  studies  of  as  reported  by  Ayele  et  al.  (2004).  Crustal  P-  and  S-wave  velocities 
below  -5  km  depth  at  FURI  range  between  5.3  km/s-6.8  km/s,  and  3.2  km/s-3.6  km/s, 
respectively  (Figure  6d).  In  addition,  beneath  FURI,  our  results  also  predict  an  ~50  km 
thick  layer  immediately  below  the  Moho  in  the  upper  mantle  that  has  P-and  S-  wave 
velocities  of  ~7.1  km/s  and  -4.3  km/s,  respectively,  which  are  anomalously  slow  (Figure 
6d).  A  similar  layer  with  P-and  S-  wave  velocities  of  -7.4  km/s  and  -4.2  km/s, 
respectively,  was  also  obtained  by  Ayele  et  al.  (2004)  (Figure  6d).  As  noted  by  Ayele  et 
al.  (2004),  this  anomalously  slow  layer  possibly  indicates  altered  lithospheric  material, 
and  supports  an  earlier  result  from  Rayleigh  wave  dispersion  of  an  approximately  100  km 
thick  lithosphere  beneath  FURI.  Station  KMBO  is  located  in  Kenya,  close  to  the  southern 
end  of  the  eastern  branch  of  the  East  African  Rift  System,  but  outside  its  edge.  Beneath 
KMBO  our  estimate  of  the  crustal  thickness  is  -38  km  (Figure  6e).  This  estimate  is 
similar  to  that  obtained  using  receiver  function  analysis  by  Dugda  et  al.  (2005)  (-41  km). 
Crustal  P-  and  S-wave  velocities  show  a  gradational  increase  with  depth  and  range 
between  5.4  km/s-8  km/s,  and  3.5  km/s-4.5  km/s,  respectively  (Figure  6e).  The  velocity 
structure  beneath  KMBO  appears  to  be  fairly  simple.  Although  it  has  relatively  high  P- 
wave  velocities  in  the  lower  crust  (Figure  6e),  it  is  otherwise  typical  of  cratonic  regions. 
Seismic  station  MBAR  is  located  between  the  western  boundary  of  the  Tanzania  craton 
and  the  western  branch  of  the  East  African  Rift  System.  Due  to  poor  correlation  of  some 
of  the  observed  phases  with  synthetics,  and  the  availability  of  only  one  event,  the 
resulting  velocity  model  from  our  study  beneath  MBAR  is  poorly  constrained. 
Nevertheless,  our  estimate  of  crustal  thickness  beneath  MBAR,  to  our  knowledge  the  first 
of  its  kind,  is  -41  km  (Figure  6f),  and  is  consistent  with  that  obtained  from  other  stations 
in  the  region.  Crustal  P-  and  S-wave  velocities  at  MBAR  range  between  5.3  km/s-7.7 
km/s,  and  3.2  km/s-3.8  km/s,  respectively  (Figure  6f).  Our  model  also  predicts  a  low  P- 
wave  velocity  (-7.2  km/s)  layer  beneath  the  crust  (Figure  6f).  Such  a  layer  promotes  the 
generation  of  the  SPL  phase  near  the  station,  which  we  observe  in  both  data  and 
synthetics  for  the  event  recorded  at  MBAR.  Therefore,  in  spite  of  the  poorly  constrained 
model  obtained  in  this  study,  we  cannot  rule  out  the  possibility  of  its  existence. 

Southern  Africa  is  another  of  the  better  studied  regions  in  Africa.  To  add  to  the 
existing  knowledge  base,  our  study  generated  P-  and  S-wave  velocity  models  beneath  the 
seismic  stations  TSUM,  LSZ,  and  LBTB  (Figures  7a-c).  Although  we  analyzed  seismic 
data  recorded  at  SUR,  due  to  the  lack  of  identifiable  phases,  we  do  not  generate  P-  and  S- 
wave  velocity  models  for  the  station.  Similar  to  most  of  the  results  in  north  and  west 
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Africa,  our  results  for  southern  Africa  are  representative  of  stable  shield  regions. 
However,  in  general,  we  obtain  slightly  higher  crustal  thicknesses  ranging  between  -42 
km  and  46  km  (Figures  7a-c).  We  also  predict  crustal  low  velocity  zones  as  discussed 
later  in  our  models  beneath  two  of  the  three  stations  in  southern  Africa.  Beneath  TSUM, 
to  our  knowledge,  no  prior  velocity  model  exists.  Thus,  the  velocity  model  obtained  from 
our  study  at  TSUM  is  the  first  of  its  kind.  We  obtained  a  crustal  thickness  beneath  TSUM 
of  -42  km  (Figure  7a).  The  velocities  of  P-  and  S-waves  in  the  crust  range  between  6.3 
km/s-7.3  km/s,  and  3.2  km/s-4  km/s,  respectively  (Figure  7a).  These  results  are  similar 
to  those  obtained  for  the  seismic  stations  in  north  and  west  Africa  and  are  therefore 
representative  of  stable  shield  regions.  We  do  not  obtain  any  anomalous  P-  and  S-wave 
velocity  zones  beneath  TSUM  (Figure  7a).  At  LSZ,  our  study  indicates  that  the  Moho  is 
located  at  a  depth  of  -43  km  (Figure  7b),  which  is  consistent  with  that  obtained  by  Midzi 
and  Ottemoller  (2001)  (-40-43  km).  The  crustal  P-wave  velocity  is  nearly  constant  (-6.2 
km/s)  between  -8  km  to  32  km  depth  (Figure  7b).  Below  -32  km  depth,  P-wave 
velocities  increase  rather  sharply  from  -6.2  km/s  to  7.8  km/s  at  the  Moho  (Figure  7b). 


pn.'ij  VSwf  |kr*ij 

Figure  7:  Obtained  P-  and  S-wave  velocity  models  for  (a)  TSUM  (b)  LSZ  and  (c)  LBTB. 
The  P-  and  S-wave  velocity  models  (broken  lines)  in  (b)  and  (c)  are  from  receiver 
function  studies  by  Midzi  and  Ottemoller  (2001). 

The  crustal  S-wave  velocities  range  between  -3.6  km/s  and  3.8  km/s  (Figure  7b).  We 
also  obtain  a  lower  crustal  low-velocity  zone  of  -5  -  8  km  thickness  in  our  velocity 
model  for  LSZ  (Figure  7b),  consistent  with  models  for  seismic  stations  elsewhere  in 
Africa  in  the  cratonic  regions.  Such  a  phenomenon  was  also  noted  by  Midzi  and 
Ottemoller  (2001).  Figure  5k  shows  the  P-  and  S-wave  velocity  models  beneath  LBTB 
obtained  in  our  study.  There  appears  to  be  a  broad  crust-mantle  transition  zone  beneath 
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LBTB  and  the  upper  bound  of  the  estimate  of  crustal  thickness  beneath  LBTB  is  ~46  km 
(Figure  7c).  Midzi  and  Ottemoller  (2001)  also  noted  the  same  and  predicted  a  crust- 
mantle  transition  zone  between  37-45  km.  The  crustal  P-  and  S-wave  velocities  beneath 
LBTB  range  between  5.8  km/s-7.5  km/s,  and  3.5  km/s-4.2  km/s,  respectively,  except  for 
a  distinct  low  P-velocity  (5.4  km/s)  zone  of  ~8  km  thickness  in  the  upper  crust  between 
10  km-20  km  depth  (Figure  7c).  Such  a  low-velocity  zone  was  also  obtained  by  Midzi 
and  Ottemoller  (2001)  at  similar  depths,  however,  the  P-wave  velocities  predicted  from 
our  study  for  this  zone  are  significantly  lower  than  those  predicted  by  Midzi  and 
Ottemdller  (2001).  Given  its  appearance  beneath  other  cratonic  seismic  stations  in  Africa, 
the  crustal  low-velocity  zone  appears  to  be  a  general  characteristic  of  the  region. 

In  summary,  Table  1  provides  the  crustal  P-  and  S-wave  velocities  and  depth  to 
the  Moho  for  each  seismic  station  in  Africa  as  obtained  from  our  study. 


Table  1:  Summary  of  crustal  P-  and  S-wave  velocities  and  depth  of  the  Moho  beneath 

Africa 


Station  name 

Station 

code 

Range  of 

r? 

(bn  s--J) 

Range  of 

Va 

(km  e-L) 

Depth  of 
MoJio 
(km) 

Tamanrasset.  Algeria 

TAM 

6.25-6.8 

3. 1-3.9 

36 

Dimbokro,  Cote  d’Ivoire 

DBIC 

6. 5-7. 3 

37-4.7 

41 

Mide  It,  Morocco 

MDT 

6. 2-6.4 

3.4-37 

38 

MTSolij;  Sei'pcgal 

MBO 

5. 6-7.2 

3.1-4. 1 

42 

Arta  Tunnel.  Dj  i  bou  ti 

AID 

4. 7-7.2 

2. 5-4. 3 

21 

Mount  Furi,  Ethiopia 

FURI 

5. 3-6. 8 

3. 2-3.6 

39 

Kilima  Mbogo,  Kenya 

KMBO 

5. 4-8.0 

3. 5-4. 5 

38 

Mbarara.  Uganda 

MB.AR 

6. 2-7.6 

2. 9-4.2 

36 

Tsumeb.  Namibia 

TSUM 

6. 3-7. 3 

3. 2-4.0 

42 

Lusaka,  Zambia 

LSZ 

6. 2-7. 8 

3. 6-3. 8 

43 

Lobatse,  Botswana 

LBTB 

5. 4-7. 5 

3. 5-4.2 

46 

4.2  China 

Based  on  our  selection  criteria,  we  analyzed  129  earthquakes  recorded  at  eleven 
permanent  broadband  seismic  stations  in  China.  The  number  of  earthquakes  recorded  at 
each  station  range  between  3  and  31.  Figure  8  shows  the  locations  of  these  seismic 
stations  and  the  geographical  distribution  of  the  earthquakes  analyzed.  The  stations 
encompass  tectonic  provinces  such  as  the  north  and  south  China  blocks,  Tibetan  plateau, 
and  Tien  Shan  Mountains.  Although  all  the  events  recorded  at  every  station  have  been 
analyzed,  work  is  ongoing  to  improve  some  of  the  modeling  results.  Therefore  in  this 
paper,  we  only  present  examples  from  those  stations  that  show  good  waveform  matches. 
Figure  9  shows  examples  of  waveform  fits  of  data  from  earthquakes  recorded  at  seismic 
stations  LSA,  WMQ,  and  BJT.  At  all  the  seismic  stations  we  observe  and  obtain  good 
matches  between  data  and  synthetics  generated  by  our  modeling  method  for  the  direct  S, 
Sp,  and  SsPmP  phases.  Except  the  direct  S  phase  which  we  observe  on  both  the 
components  at  all  the  stations,  we  note  that  the  Sp  phase  is  prominent  on  the  vertical  and 
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radial  components  at  LSA  and  WMQ  but  not  on  the  radial  component  at  BJT.  Similarly, 
we  observe  the  SsPmP  phase  on  both  components  at  LSA  and  WMQ  but  not  on  the 
vertical  component  at  BJT.  However,  the  SPL  phase  is  only  noted  and  well  matched  on 
the  radial  component  at  LSA.  We  confirm  the  presence  of  the  SPL  phase  by  analyzing  the 
particle  motion  within  the  corresponding  time  window  which  turns  out  to  be  prograde 
elliptical. 

240°  0°  120° 


240° 


0° 


120° 


Figure  8:  Map  showing  earthquakes  analyzed  in  this  study  (blue  stars)  and  the  permanent 
broadband  seismic  stations  in  China  (red  triangles)  that  recorded  them.  The  respective 
station  codes  are  shown  adjacent  to  location  of  each  station. 


Figure  9:  Vertical  and  radial  component  seismograms  for  example  events  recorded  at 
LSA,  WMQ,  and  BJT  showing  the  observed  (solid  line)  and  synthetic  (dashed  line) 
waveforms.  The  correlated  waveforms  are  indicated  on  the  panels. 


Based  on  waveform  fits  obtained  for  each  source-receiver  pair  we  generate 
velocity  models  for  each  seismic  station.  However,  the  models  generated  for  the  same 
station  using  waveform  fits  from  different  earthquakes  recorded  at  that  station,  although 
similar,  are  not  identical.  To  analyze  which  velocity  model  is  more  reliable  we  calculate 
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the  PPD  and  correlation  matrices  for  the  modeling  results  of  each  source-receiver  pair. 
Here,  we  show  examples  of  these  statistical  calculations  from  an  event  recorded  at  LSA 
(Figure  10a  and  b),  and  WMQ  (Figure  10c  and  d). 


LSA  parameter  torretsdoiva 


WMQ  parameter  correlations 


5  10  15  20  25  30  35  40 

parameter 


(c)  (d) 

Figure  10:  Posterior  Probability  Distribution  (PPD)  for  example  events  recorded  at  (a) 
LSA  and  (c)  WMQ.  More  peaked  distributions  indicate  more  uniqueness  among  different 
models  and  fewer  trade-offs  among  model  parameters.  Model  parameter  correlation 
matrices  for  the  same  event  at  (b)  LSA  and  (d)  WMQ  are  also  shown.  Each  small  square 
represents  a  model  parameter  (Yp,  Vs,  Thickness  of  Layer,  and  Density)  on  both  axes. 
The  correlations  range  between  -1  and  1.  Sparse  colored  squares  off-diagonal  indicate 
better  constraints  and  greater  confidence  (less  trade-offs)  in  those  parts  of  the  models. 
Note  that  at  LSA  where  we  observe  and  match  the  SPL  phase  there  are  less  colored 
squares  in  the  correlation  matrix  in  the  lower  crust-upper  mantle  (b)  compared  to  that  at 
WMQ  (d),  suggesting  that  SPL  improves  constraints  in  those  parts  of  the  models. 


Finally,  in  Figure  11  we  show  the  velocity  models  for  each  seismic  station  in 
China  that  we  obtain  from  our  modeling  exercises  so  far.  Wherever  available  we  also 
show  the  models  from  an  earlier  study  for  comparison.  Our  models  are  consistent  with 
regional  tectonics  and  models  obtained  from  earlier  studies  using  receiver  functions  and 
seismic  tomography.  The  crustal  thicknesses  beneath  stations  in  north  China  block  range 
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between  ~32  and  42  km  with  the  crust  in  central  part  of  the  block  near  station  BJT  being 
thicker.  However,  within  north  China  block  towards  its  western  edge,  the  crust  appears 
anomalously  thick  (—55  km)  beneath  station  LZH.  This  region  also  coincides  with  the 
border  of  north  China  block  and  Tibetan  plateau.  The  south  China  block  consists  of 
widely  varying  crustal  thicknesses:  ~33  km  beneath  station  ENH  in  the  northern  part  of 
the  block  and  ~52  km  beneath  station  KMI  in  the  southern  part,  for  example.  The 
southern  part  of  south  China  block  near  station  KMI  also  coincides  with  its  border  with 
the  Tibetan  plateau,  where  the  Moho  is  significantly  deeper  than  elsewhere  in  China, 
implying  a  deep  crustal  root.  Station  LSA,  in  Lhasa,  Tibet,  has  a  crustal  thickness  of  ~53 
km.  In  Tien  Shan  Mountains  in  northwestern  China,  the  crust  is  thicker  than  average  (~42 
km)  but  not  as  thick  as  that  beneath  the  Tibetan  plateau. 
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Figure  11:  Obtained  P-  and  S-wave  velocity  models  (solid  lines)  for  the  eleven 
permanent  broadband  seismic  stations  in  China.  Station  codes  are  indicated  in  each  panel. 
The  models  (broken  lines)  in  LZH,  BJT,  HIA,  MDJ,  WMQ,  and  KMI  are  from  receiver 
function  studies  by  Mangino  et  al.  (1999). 

In  summary,  Table  2  provides  the  crustal  P-  and  S-wave  velocities  and  depth  to 
the  Moho  for  each  seismic  station  in  China  as  obtained  from  our  study. 
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Table  2:  Summary  of  crustal  P-  and  S-wave  velocities  and  depth  of  the  Moho  beneath 

China 


Station 

Range  of  Vp  (km/s) 

Range  of  Vs  (km/s) 

Depth  of  Moho 
(km) 

LZH 

4. 3-8. 5 

2. 3-4. 7 

55 

BJT 

3. 8-7. 3 

2.2-4. 3 

42 

HIA 

6. 0-8.4 

3.2-4. 8 

40 

MDJ 

5.4-8. 0 

3.2-4. 3 

39 

WMQ 

4. 0-8.2 

2. 5-4.5 

42 

SSE 

4. 1-7.5 

2. 1-4.3 

39  (?) 

KMI 

6. 1-8.3 

3. 3-4. 6 

52 

LSA 

6. 0-7. 7 

3. 1-4.2 

53 

ENH 

4.0-7. 1 

2. 3-4.0 

33 

XAN 

3. 9-7. 5 

2. 3-3. 8 

32 

QIZ 

6. 0-8. 7 

3.4-4. 6 

52 

(?)  The  depth  of  IV 

[oho  computations  are  preliminary  and  requires  assessment  of 

uncertainties 

4.3  Canada 

We  apply  our  modeling  method  to  data  from  137  earthquakes  recorded 
teleseismically  at  11  permanent  broadband  seismic  stations  spanning  Canada  during 
1976-2005  (Figure  12).  We  show  examples  of  waveform  fitting  for  some  events,  and  also 
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Figure  12:  Map  showing  earthquakes  analyzed  in  this  study  (blue  stars)  and  the 
permanent  broadband  seismic  stations  in  Canada  (red  triangles)  that  recorded  them.  The 
respective  station  codes  are  shown  adjacent  to  location  of  each  station. 
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describe  interpretations  of  the  uncertainty  computations.  The  number  of  earthquakes 
recorded  at  each  station  range  between  3  and  26.  The  stations  encompass  tectonic 
provinces  such  as  Cordilleran  orogen,  western  plains  and  Slave  province  in  western 
Canada,  Grenville  province  and  Appalachian  orogen  in  eastern  Canada,  and  the  Canadian 
Arctic.  Figure  13a  shows  examples  of  waveform  fits  of  data  from  earthquakes  recorded  at 
seismic  stations  INK,  LLLB,  and  GAC.  At  all  the  seismic  stations  we  observe  and  obtain 
good  matches  between  data  and  synthetics  of  the  direct  S  and  SsPmP  phases.  We  observe 
and  match  the  Sp  phase  on  both  the  components  at  INK  but  only  on  the  vertical 
component  at  LLLB,  and  radial  component  at  GAC.  The  SPL  phase  appears  abundant  at 
the  Canadian  seismic  stations  and  we  observe  and  obtain  good  matches  on  the  vertical 
component  at  INK,  both  components  at  LLLB,  and  radial  component  at  GAC.  We 
confirm  our  observation  of  the  SPL  phase  by  observing  prograde  elliptical  particle 
motion  diagrams  (Figure  13b).  We  also  calculate  the  posterior  probability  distribution 
(PPD)  and  correlation  matrix  for  each  source-receiver  pair,  examples  of  which  we  show 
in  Figure  13c.  In  Figure  14  we  present  the  final  velocity  models  for  each  seismic  station 
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Figure  13:  (a)  Vertical  and  radial  component  seismograms  for  example  events  recorded 
at  INK,  LLLB,  and  GAC  showing  the  observed  (solid  line)  and  synthetic  (dashed  line) 
waveforms.  The  correlated  waveforms  are  indicated  on  the  panels,  (b)  Particle  motion 
diagrams  for  a  time  window  of  8  seconds  around  the  SPL  phase  on  the  data  and 
synthetics  for  an  event  recorded  at  INK,  and  10  seconds  at  LLLB  and  GAC  showing 
prograde  elliptical  motion  diagnostic  of  the  SPL  phase.  The  red  portions  of  the  diagrams 
indicate  beginning  of  the  motion. 
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Figure  13:  (c)  Marginal  Posterior  Probability  Distributions  (PPD)  for  example  events 
recorded  at  seismic  stations  DRLN  and  GAC.  More  peaked  distributions  indicate  more 
uniqueness  among  different  models  and  fewer  trade-offs  among  model  parameters. 
Model  parameter  correlation  matrices  for  the  same  events  at  DRLN  and  GAC.  Each  small 
square  represents  a  model  parameter  on  both  axes.  The  correlations  range  between  -1  and 
1 .  Sparse  colored  squares  off-diagonal  indicate  better  constraints  and  lesser  trade-offs  in 
those  parts  of  the  models. 


in  Canada.  Wherever  available,  we  also  show  the  earlier  velocity  models  (Cassidy,  1995; 
Darbyshire,  2003)  for  comparison.  Our  models  beneath  the  Canadian  seismic  stations  are 
consistent  with  earlier  studies  (e.g.  Ramesh  et  al.,  2002)  and  regional  tectonics  (Figure 
14).  Crustal  thicknesses  beneath  stations  in  the  northern  Cordilleran  orogen,  western 
plains,  and  Slave  province  range  between  35  and  37  km.  The  Moho  appears  to  be  slightly 
shallower  (30-36  km)  beneath  stations  in  the  southern  Cordilleran  orogen.  In  eastern 
Canada,  the  crust  beneath  stations  of  the  Grenville  province  and  Appalachian  orogen  is 
generally  thick  (—44  km)  with  the  exception  of  that  beneath  station  DRLN,  in  the 
northeast,  where  it  is  ~33  km.  Moho  depths  beneath  stations  in  the  Canadian  Arctic  range 
between  ~30  and  41  km.  We  also  observe  low- velocity  zones  in  the  crust  and  uppermost 
mantle  at  some  stations,  however  the  constraints  on  them  are  not  strong.  Taken  together, 
the  results  provide  a  comprehensive  snapshot  of  the  velocity  structure  beneath  Canada. 


26 


Figure  14:  P-  and  S-wave  velocity  models  (solid  lines)  for  seismic  stations  in  Canada. 
The  P-  and  S-wave  velocity  models  in  FRB,  MBC,  and  RES  are  from  receiver  function 
studies  by  Darbyshire  (2003),  and  in  GAC  and  INK  are  those  from  Cassidy  (1995). 
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In  summary,  Table  3  provides  the  crustal  P-  and  S-wave  velocities,  depth  to  the 
Moho,  and  average  crustal  Poisson’s  ratio  for  each  seismic  station  in  Canada  as  obtained 
from  our  study.  (P.S.  The  Poisson’s  ratio  computations  are  preliminary) 

Table  3:  Summary  of  crustal  P-  and  S-wave  velocities,  depth  of  the  Moho,  and  average 

crustal  Poisson’s  ratio  beneath  Canada 


Station 

Range  of  Vp 
(km/s) 

Range  of  Vs 
(km/s) 

Depth  of 
Moho  (km) 

Av.  Crustal 
Poisson’s 
Ratio 

DRLN 

6. 1-8.1 

3. 3-4. 3 

33 

0.24 

FRB 

6. 0-7. 8 

3. 3-4. 8 

41 

0.25 

GAC 

6. 5-7. 6 

3.4-4. 6 

44 

0.29 

INK 

5. 6-8. 7 

3.2-4.4 

35 

0.21 

LLLB 

62-1.5 

3.2-4. 3 

36 

0.20 

MBC 

42-1.2 

2.4-4. 0 

30 

0.22 

PMB 

5. 8-7. 9 

3. 4-4.6 

30 

0.24 

RES 

3. 6-8. 5 

2.2-4. 6 

38 

0.26 

WHY 

6.5-1.5 

3. 3-4.4 

37  (?) 

0.24 

YKW3 

6.0-8. 1 

3.4-4.4 

35 

0.24 

SCHQ 

6.6-12 

3.2-4.4 

31  (?) 

0.27 

(?)  The  depth  of  Moho  computations  are  preliminary  and  requires  assessment  of 
uncertainties 


5.  CONCLUSIONS 

In  this  report,  we  discussed  a  waveform  fitting  technique  that  relies  on  a 
parallelized  reflectivity  method  to  compute  synthetic  seismograms  and  implements  a 
global  optimization  algorithm  using  Very  Fast  Simulated  Annealing  (VFSA).  We  also 
demonstrated  the  application  of  the  method  to  determine  one-dimensional,  azimuthally- 
dependent,  crust  and  upper  mantle  P-  and  S-wave  velocity  structure  beneath  broadband 
seismic  stations  across  the  continent  of  Africa,  China,  and  Canada.  Our  technique  is 
complementary  to  receiver  function  methods,  has  many  of  its  advantages,  and  adds  to 
them  in  providing  regional  crust  and  upper  mantle  structure  appropriate  for  locating 
small,  regional  events.  We  are  also  able  to  compute  synthetic  seismograms  that  contain 
all  the  possible  phases  for  a  prescribed  source-receiver  path,  and  obtain  direct  estimates 
of  the  P-  and  S-wave  velocities  beneath  seismic  stations.  Statistical  tools  incorporated  in 
the  technique  allow  us  to  assess  uncertainties  associated  with  our  models  and  estimate 
tradeoffs  between  model  parameters  in  different  layers.  The  use  of  the  SPL  phase  as 
shown  in  the  study,  enhances  our  constraints  for  lower  crust  and  upper  mantle  structure 
beneath  the  seismic  stations. 

Advantages  of  the  waveform  modeling  method  described  here  include  the  ability 
to  model  simultaneously  all  phases  that  might  be  present  in  the  observed  waveform, 
provide  independent  estimates  of  both  P-  and  S-wave  velocities,  assess  uncertainties  in 
model  parameters,  find  a  range  of  acceptable  models  that  explain  the  data,  and  minimize 
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dependence  of  the  final  results  on  the  initial  model.  Furthermore,  a  full  waveform  fitting 
procedure  iterated  by  VFSA  reduces  the  opportunities  for  imposing  bias  by  automating 
the  choice  of  model  perturbations  and  by  eliminating  the  need  for  a  user  to  pick 
functionals  of  the  seismogram  (e.g.,  arrival  time,  amplitude)  precisely  or  accurately,  or  to 
identify  phases  correctly  and  construct  associated  raypaths.  The  procedure  is  therefore 
relatively  free  of  user  bias  and  robust  compared  to  methods  in  which  one  measures 
functionals  in  pre-processing  steps.  In  our  method,  we  avoid  the  risk  of  incorrectly 
identifying  phases  or  incorrectly  measuring  arrival  times  due  to  wave  interference  or  pure 
blunders.  All  phases  possible  for  each  iteration’s  model  will  be  included  in  its 
corresponding  synthetic  seismograms  (including  appropriate  interference  and  distortion), 
they  will  automatically  be  evaluated  against  the  observed  data  without  risk  that  mis- 
identification  will  map  mis-fits  to  incorrect  raypaths,  for  example. 

Applied  to  large-magnitude,  deep-focus  earthquakes  recorded  teleseismically  in 
Africa,  China,  and  Canada  our  method  successfully  produced  crust  and  upper  mantle 
(wherever  SPL  was  observed)  P-  and  S-wave  velocity  models,  that  are  consistent  with 
earlier  models,  in  the  sense  that  they  fall  within  the  associated  uncertainties  we  found 
with  the  products  of  multiple  VFSA  runs.  For  some  seismic  stations  in  Africa,  our  study 
provided  such  velocity  models  that  are  the  first  of  their  kind.  Our  models  were  also 
consistent  with  the  broad  regional  tectonics  of  the  respective  regions.  While  the  technique 
described  here  provided  layered,  one-dimensional  models,  a  dataset  that  includes  a 
broader  azimuthal  distribution  of  earthquakes  for  each  station  would  allow  this  source- 
receiver-based  technique  to  produce  better  azimuthally-dependent  models,  and  thus  a 
more  detailed  view  of  the  Earth  structure. 


29 


30 


REFERENCES 


Ammon,  C.  J.,  Randall,  G.  E.  &  Zandt,  G.,  1990.  On  the  nonuniqueness  of  receiver 
function  inversions,  J.  Geophys.  Res.,  95,  15,303  -  15,318. 

Ayele,  A.,  Stuart,  G.  &  Kendall,  J.  -M.,  2004.  Insights  into  rifting  from  shear  wave 
splitting  and  receiver  functions:  an  example  from  Ethiopia,  Geophys.  J.  Int., 
157,354-362. 

Baag,  C.  E.  &  Langston,  C.  A.,  1985.  Shear-coupled  PL,  Geophys.  J.  Int.,  80,  363  -  385. 

Baag,  C.  E.  &  Langston,  C.  A.,  1986.  Diffracted  Sp  Generated  Under  the  Australian 
Shield,  J.  Geophys.  Res.,  91,  9,507  -  9,516. 

Bock,  G.,  1988.  Sp  phases  from  the  Australian  upper  mantle,  Geophys.  J.  Int.,  94,  73  - 
81. 

Bock,  G.,  1991.  Long-Period  S  To  P  Converted  Waves  And  the  Onset  of  Partial  Melting 
Beneath  Oahu,  Hawaii,  Geophys.  Res.  Lett.,  18,  869  -  872. 

Bock,  G.  &  Kind,  R.,  1991.  A  global  study  of  S-to-P  and  P-to-S  conversions  from  the 
upper  mantle  transition  zone,  Geophys.  J.  Int.,  107,  1 17  -  129. 

Cassidy,  J.  F.,  1995.  A  comparison  of  the  receiver  structure  beneath  stations  of  the 
Canadian  National  Seismic  Network,  Can,  J.  Earth  Sc.,  32,  938-951. 

Chander,  R.,  Alsop,  L.  E.  &  Oliver,  J.,  1968.  On  the  synthesis  of  shear  coupled  PI  waves, 
Bull.  Seismol.  Soc.  Am.,  58,  1,849  -  1,877. 

Darbyshire,  F.  A.,  2003.  Crustal  structure  across  the  Canadian  High  Arctic  region  from 
teleseismic  receiver  function  analysis,  Geophys.  J.  Int.,  152,  372-391. 

Dugda,  M.  T.,  Nyblade,  A.  A.,  Julia,  J.,  Langston,  C.  A.,  Ammon,  C.  J.  &  Simiyu,  S., 
2005.  Crustal  structure  in  Ethiopia  and  Kenya  from  receiver  function  analysis: 
Implications  for  rift  development  in  eastern  Africa,  J.  Geophys.  Res.,  110, 
B01303,  doi:10.1029/2004JB003065. 

Dugda,  M.  T.  &  Nyblade,  A.  A.,  2006.  New  constraints  on  crustal  structure  in  eastern 
Afar  from  the  analysis  of  receiver  functions  and  surface  wave  dispersion  in 
Djibouti,  in  Yirgu,  G.,Ebinger,  C.  J.  &  Maguire,  P.  K.  H.  (eds.),  The  Structure 
and  Evolution  of  the  East  African  Rift  System  in  the  Afar  Volcanic  Province, 
Geol.  Soc.  London,  Special  Publications,  259,  241  -253. 

Dziewonski,  A.  M.  &  Anderson,  D.  L.,  1981.  Preliminary  reference  Earth  model,  Phys. 
Earth  planet.  Int.,  25,  297  -  356. 

Frazer,  L.  N.,  1977.  Synthesis  of  shear  coupled  PL,  Ph.D  Thesis,  Princeton  University, 
Princeton,  New  Jersey,  U.S.A.,  54  pp. 

Gangopadhyay,  A.,  Pulliam,  J.,  Sen,  M.  K.,  (2007).  Waveform  Modeling  of  Teleseismic 
S,  Sp,  SsPmP,  and  Shear-Coupled  PL  waves  for  Crust  and  Upper  Mantle  Velocity 
Structure  Beneath  Africa.  Geophys.  J.  Int.,  170,  1,210  -  1,226, 

doi:  10.1 1 1 1/j.  1365-246X.2007.03470.X. 

Gropp,  W.  &  Lusk,  E.,  1995.  Dynamic  process  management  in  an  MPI  setting,  7th  IEEE 
Symposium  on  Parallel  and  Distributeed  Processing,  p.  530. 

Hazier,  S.  E.,  Sheehan,  A.  F.,  McNamara,  D.  E.  &  Walter,  W.  R.,  2001.  One-dimensional 
Shear  Velocity  Structure  of  Northern  Africa  from  Rayleigh  Wave  Group  Velocity 
Dispersion,  Pageoph.,  158,  1,475  -  1,493. 

Ingber,  L.,  1989.  Very  fast  simulated  reannealing,  Mathl.  Comput.  Modeling,  12,  967  - 
993. 


31 


Jordan,  T.  H.  &  Frazer,  L.  N.,  1975.  Crustal  and  upper  mantle  structure  from  Sp  phases, 
J.  Geophys.  Res.,  80,  1,504-  1,518. 

Julia,  J.,  Ammon,  C.  J.,  Herrmann,  R.  B.  &  Correig,  A.  M.,  2000.  Joint  Inversion  of 
Receiver  Function  and  Surface  Wave  Dispersion  Observations,  Geophys.  J.  Int., 
143,99-112. 

Kennett,  B.  L.  N.,  1983.  Seismic  wave  propagation  in  stratified  media,  Cambridge 
University  Press,  Cambridge,  338  pp. 

Langston,  C.  A.,  1996.  The  SsPmP  Phase  in  Regional  Wave  Propagation,  Bull.  Seismol. 
Soc.  Am.,  86,  133  -  143. 

Mallick,  S.  &  Frazer,  L.  N.,  1987.  Practical  aspects  of  reflectivity  modeling,  Geophysics, 
52,1,355-  1,364. 

Mangino,  S.,  Priestley,  K.  &  Ebel,  J.,  1999.  The  Receiver  structure  beneath  the  China 
Digital  Seismograph  stations,  Bull.  Seis.  Soc.  Am.,  89,  1,053-1,076. 

Marone,  F.,  van  der  Meijde,  M.,  van  der  Lee,  S.  &  Giardini,  D.,  2003.  Joint  inversion  of 
local,  regional  and  teleseismic  data  for  crustal  thickness  in  the  Eurasia-Affica 
plate  boundaryregion,  Geophys.  J.  Int.,  154,  499  -  514. 

Midzi,  V.  &  Ottemoller,  L.,  2001.  Receiver  function  structure  beneath  three  southern 
Africa  seismic  broadband  stations,  Tectonophysics,  339,  443  -  454. 

Oliver,  J.,  1961.  On  the  long  period  character  of  shear  waves,  Bull.  Seismol.  Soc.  Am.,  51, 
1-12. 

Owens,  T.  J.,  Zandt,  G.  &  Taylor,  S.  R.,  1984.  Seismic  evidence  for  an  ancient  rift 
beneath  the  Cumberland  Plateau,  Tennessee:  A  detailed  analysis  of  broadband 
teleseismic  P  waveforms,  J.  Geophys.  Res.,  89,  7,783  -  7,795. 

Owens,  T.  J.  &  Zandt,  G.,  1997.  Implications  of  crustal  property  variations  for  models  of 
Tibetan  plateau  evolution,  Nature,  387,  37-43. 

Pasyanos,  M.  E.  &  Walter,  W.  R.,  2002.  Crust  and  upper-mantle  structure  of  North 
Africa,  Europe  and  the  Middle  east  from  inversion  of  surface  waves,  Geophys.  J. 
Int.,  149,463-481. 

Pasyanos,  M.  E.,  Walter,  W.  R.,  Flanagan,  M.  P.,  Goldstein,  P.  &  Bhattacharyya,  J., 
2004.  Building  and  Testing  an  a  priori  Geophysical  Model  for  Western  Eurasia 
and  North  Africa,  Pageoph.,  161,  235  -  281. 

Pasyanos,  M.  E.,  2005.  A  variable  resolution  surface  wave  dispersion  study  of  Eurasia, 
North  Africa,  and  surrounding  regions,  J.  Geophys.  Res.,  110,  B12301, 
doi:  10. 1029/2005 JB003749. 

Pulliam,  J.,  Sen,  M.  K.,  Frohlich,  C.  &  Grand,  S.,  2002.  Waveform  Modeling  of  the  Crust 
and  Upper  Mantle  Using  S,  Sp,  SsPmP  and  Shear-Coupled  PL  waves,  Proc.  24th 
Seis.  Res.  Rev.  -  Nuclear  Explosion  Monitoring:  Innovation  and  Integration,  144 
-  153. 

Pulliam,  J.  &  Sen,  M.,  2004.  Waveform  Modeling  of  the  Crust  and  Upper  Mantle  using 
S,  Sp,  SsPmP,  and  Shear-Coupled  PL  waves  for  improved  event  location,  focal 
depth  determination,  and  uncertainty  estimation,  Proc.  26th  Seis.  Res.  Rev.  - 
Trends  in  Nuclear  Explosion  Monitoring,  142  -  152. 

Pulliam,  J.  &  Sen,  M.  K.,  2005.  Assessing  Uncertainties  in  Waveform  Modeling  of  the 
Crust  and  Upper  Mantle,  Proc.  27th  Seis.  Res.  Rev.  -  Ground-Based  Nuclear 
Explosion  Monitoring  Technologies,  152  -  160. 


32 


Ramesh,  D.  S.,  Kind,  R.  &  Yuan,  X.,  2002.  Receiver  function  analysis  of  the  North 
American  crust  and  upper  mantle,  Geophys.  J.  Int.,  150,  91  -  108. 

Sandvol,  E.,  Seber,  D.,  Calvert,  A.  &  Barazangi,  M.,  1998.  Grid  search  modeling  of 
receiver  functions:  Implications  for  crustal  structure  in  the  Middle  East  and  North 
Africa,  J.  Geophys.  Res.,  103,  26,899  -  26,917. 

Sen,  M.  K.  &  Stoffa,  P.  L.,  1991.  Non-linear  one-dimensional  seismic  waveform 
inversion  using  simulated  annealing,  Geophysics ,  56,  1624  -  1638. 

Sen,  M.  K.  &  Stoffa,  P.  L.,  1995.  Global  Optimization  Methods  in  Geophysical 
Inversion,  Elsevier  Science  Publishing  Company,  The  Netherlands,  281  pp. 

Sen,  M.  K.  &  Stoffa,  P.  L.,  1996.  Bayesian  inference,  Gibbs’  sampler  and  uncertainty 
estimation  in  geophysical  inversion,  Geophysical  Prospecting,  44,  313  -  350. 

Swenson,  J.  L.,  Beck,  S.  L.  &  Zandt,  G.,  1999.  Regional  distance  shear-coupled  PL 
propagation  within  the  northern  Altiplano,  central  Andes,  Geophys.  J.  Int.,  139, 
743  -  753. 

Tarantola,  A.,  1994.  Inverse  Problem  Theory:  Methods  for  Data  Fitting  and  Model 
Parameter  Estimation,  Elsevier,  Amsterdam,  The  Netherlands,  600  pp. 

Tiberi,  C.,  Ebinger,  C.,  Ballu,  V.,  Stuart,  G.  &  Oluma  B.,  2005.  Inverse  models  of  gravity 
data  from  the  Red  Sea  -  Aden  -  East  African  rifts  triple  junction  zone,  Geophys. 
J.  Int.,  163,  775  -  787. 

Wessel,  P.  &  Smith,  W.  H.  F.,  1991.  Free  software  helps  map  and  display  data,  EOS 
Trans.  AGU,  72,  441. 

Zandt,  G.  &  Randall,  G.  E.,  1985.  Observations  of  Shear-Coupled  P  Waves,  Geophys. 
Res.  Lett.,  12,  565  -  568. 

Zandt,  G.,  Beck,  S.  L.,  Ruppert,  S.  R.,  Ammon,  C.  J.,  Rock,  D.,  Minaya,  E.,  Wallace,  T. 
C.  &  Silver,  P.  G.,  1996.  Anomalous  Crust  of  the  Bolivian  Altiplano,  Central 
Andes:  Constraints  from  Broadband  Regional  Seismic  Waveforms,  Geophys.  Res. 
Lett.,  23,  1,159  -  1,162. 

Zhang,  J.  &  Langston,  C.  A.,  1996.  Array  Observations  of  the  Shear-Coupled  PL  Wave, 
Bull.  Seismol.  Soc.  Am.,  86,  538  -  543. 

Zhao,  L.  -S.  &  Frohlich,  C.,  1996.  Teleseismic  body-waveforms  and  receiver  structures 
beneath  seismic  stations,  Geophys.  J.  Int.,  146,  525  -  540. 

Zhao,  L.  -S.,  Sen,  M.  K.,  Stoffa,  P.  L.  &  Frohlich,  C.,  1996.  Application  of  very  fast 
simulated  annealing  to  the  determination  of  the  crustal  structure  beneath  Tibet, 
Geophys.  J.  Int.,  125,  355  -  370. 


33 


34 


IRIS 

GMT 

VFSA 

PREM 

SPL 

MPI 

SA 

PPD 

CMT 


List  of  Symbols,  Abbreviations,  and  Acronyms 

Incorporated  Research  Institutions  for  Seismology 

Generic  Mapping  Tools 

Very  Fast  Simulated  Annealing 

Preliminary  Reference  Earth  Model 

Shear-Coupled  PL 

Message  Passing  Interface 

Simulated  Annealing 

Posterior  Probability  Density 

Centroid  Moment  Tensor 


35 


